﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace DynaPrecision.LinearSystem
{
	/// <summary>
	/// LU Decomposition-based solver, a Gauss Elimination method
	/// </summary>
	/// <typeparam name="TData"></typeparam>
	public class DoolittleDecompositionSolver: ILinearSystemSolver, IBandArraySolver
	{

		#region ILinearSystemSolver Members
		/// <summary>
		/// Solve assuming a general array; the decomposed array will be stored in the memory. When RAM is limited, considering first Domposing then call the decomposed array solver
		/// </summary>
		/// <typeparam name="TInput"></typeparam>
		/// <typeparam name="TLoad"></typeparam>
		/// <typeparam name="TResult"></typeparam>
		/// <param name="input"></param>
		/// <param name="load"></param>
		public void SolveFor<TData>(IMultiDimensionalArray<TData> input, IMultiDimensionalArray<TData> load, IMultiDimensionalArray<TData> result)
		{
			var tempStorage = new MultiDimensionalArray<TData>(input.GetLength(0), input.GetLength(1));
			DoolittleDecomposition.Decompose(input, tempStorage);
			for (int i = 0; i < tempStorage.GetLength(0); i++)
			{
				result[i, 0] = load[i, 0];
				for (int k = 0; k < i; k++)
				{
					result[i, 0] = result[i, 0] - result[k, 0] * (dynamic)tempStorage[i, k];
				}
			}
			for (int i = tempStorage.GetLength(0) - 1; i >= 0; i--)
			{
				for (int k = tempStorage.GetLength(0) - 1; k > i; k--)
				{
					result[i, 0] = result[i, 0] - result[k, 0] * (dynamic)tempStorage[i, k];
				}
				result[i, 0] = result[i, 0] / (dynamic)tempStorage[i, i];
			}
			return;
		}

		#endregion

		#region IBandArraySolver<TInput> Members

		public void SolveFor<TData>(IMultiDimensionalArray<TData> input, IMultiDimensionalArray<TData> load, IMultiDimensionalArray<TData> result, int halfBandWidth)
		{
			BandArray<TData> tempStorage = new BandArray<TData>(input.GetLength(0), halfBandWidth);
			DoolittleDecomposition.Decompose(input, tempStorage, halfBandWidth);
			for (int i = 0; i < tempStorage.GetLength(0); i++)
			{
				result[i, 0] = load[i, 0];
				int k = 0 > i - halfBandWidth ? 0 : i - halfBandWidth;
				for (; k < i; k++)
				{
					result[i, 0] = result[i, 0] - result[k, 0] * (dynamic)tempStorage[i, k];
				}
			}
			for (int i = tempStorage.GetLength(0) - 1; i >= 0; i--)
			{
				int k = tempStorage.GetLength(0) - 1 > i + halfBandWidth ? i + halfBandWidth : tempStorage.GetLength(0) - 1;
				for (; k > i; k--)
				{
					result[i, 0] = result[i, 0] - result[k, 0] * (dynamic)tempStorage[i, k];
				}
				result[i, 0] = result[i, 0] / (dynamic)tempStorage[i, i];
			}
			return;
		}
		#endregion
	}
}